'''
Created on 2010-11-9

@author: yun_hua
'''

import os
import unittest
import numpy as N
import pylab as pl

from numheat.tubreactor import TubularReactor
from numheat.geo import Grid
from numheat.paramfactory import ParamFactory
from numheat.solver import JacobiSolver, GaussSeidelSolver

class TubReactorTest(unittest.TestCase):
    """
    Test for the tubular reactor.
    """
    def testcstr(self):
        param = ParamFactory()
        reactor = TubularReactor(Grid.rectangle(10., 0.5, 40, 10), 'Upwind', param)
        self.assertAlmostEqual(reactor.get_deltax(), 0.25)
        self.assertAlmostEqual(reactor.get_deltay(), 0.05)
        
        # No mass or heat diffusivities are given.
        self.assertEqual(reactor.getmarray().shape, (10 + 2, 40 + 2))
        self.assertEqual(reactor.gettarray().shape, (10 + 2, 40 + 2))
        #self.assertAlmostEqual(reactor.getaE(5, 4, False), 0.175 * 1.0 * 0.05 / 0.25)
        #self.assertAlmostEqual(reactor.getaW(5, 4, False), 0.175 * 1.0 * 0.05 / 0.25)
        
    def testrun(self):
        radius = 0.1
        length = 1.
        nx = 100
        ny = 40
        param = ParamFactory()
        reactor = TubularReactor(Grid.rectangle(length, radius, nx, ny), 'Upwind', param)
        s = GaussSeidelSolver(reactor, delta=1e-6, maxiter=100, verbosity=True)
        s.run()
        
        #print s.getsol()
        #print reactor.gettarray()
        xi = reactor.getgrid()[0, :, 0]
        yi = reactor.getgrid()[:, 0, 1]
        zimass = s.getsol()[0]
        zitemp = s.getsol()[1]
        conv = s.getsol()[2]
        con = s.getsol()[3]
        self.assertEqual(len(xi), len(conv))
        self.assertTrue(conv[0] < 1.)
        
        pl.matplotlib.pyplot.autumn()
        save_path = "data"
        N.save(os.path.join(save_path, "Mass%03d") % 1, zimass)
        N.save(os.path.join(save_path, "Heat%03d") % 1, zitemp)     
        
        #mass
        pl.contourf(xi, yi, zimass, 100)
        cont = pl.contour(xi, yi, zimass, 10, linewidths=0.5, colors='blue')
        pl.clabel(cont, inline=1, fontsize=8)
        pl.axis([0, length, 0, radius])
        pl.title("Mass diffusion in FT reactor.")
        pl.xlabel("Length (m)")
        pl.ylabel("Radius (m)")
        cbar = pl.colorbar(cont)
        cbar.add_lines(cont)
        save_path = "data"
        pl.savefig(os.path.join(save_path, "Mcontour%03d.png") % 1)
        pl.clf()
        
        # temperature
        pl.contourf(xi, yi, zitemp, 10)
        cont = pl.contour(xi, yi, zitemp, 10, linewidths=0.5, colors='k')
        pl.clabel(cont, inline=1, fontsize=8)
        pl.axis([0, length, 0, radius])
        pl.title("Heat diffusion in FT reactor.")
        pl.xlabel("Length (m)")
        pl.ylabel("Radius (m)")
        cbar = pl.colorbar(cont)
        cbar.add_lines(cont)
        pl.savefig(os.path.join(save_path, "Tcontour%03d.png") % 1)
        pl.clf()
        
        # Conversion profile
        pl.title("Average conversion")
        pl.plot(xi, conv, '-', label="X_A", linewidth=1)
        pl.xlabel('len')
        pl.ylabel('conversion')
        pl.legend()
        pl.savefig(os.path.join(save_path, "Conv%03d.png") % 1)
        pl.clf()
        
        # Concentration profile
        pl.title("Mean concentration(mol/m3)")
        pl.plot(xi, con, '-', label="C_A", linewidth=1)
        pl.xlabel('len')
        pl.ylabel('concentration')
        pl.legend()
        pl.savefig(os.path.join(save_path, "ConA%03d.png") % 1)
        pl.clf()